subdiagnosis <- readr::read_tsv(
file.path("..", "..", "..", "data", "current", params$scpca_project_id, "single_cell_metadata.tsv"),
show_col_types = FALSE
) |>
dplyr::filter(scpca_sample_id == params$sample_id) |>
dplyr::pull(subdiagnosis)This notebook explores using CopyKAT
to estimate tumor and normal cells in SCPCS000194 from SCPCP000006. This
sample has a(n) Anaplastic subdiagnosis.
CopyKAT was run using the copyKAT.R script
with and without a normal reference, using euclidean to calculate
distance in copyKAT.
These results are read into this notebook and used to:
CopyKAT to annotate tumor cells.CopyKAT to cell type
annotations using label transfer and the fetal (kidney) references.# The base path for the OpenScPCA repository, found by its (hidden) .git directory
repository_base <- rprojroot::find_root(rprojroot::is_git_root)
# The current data directory, found within the repository base directory
data_dir <- file.path(repository_base, "data", "current", params$scpca_project_id)
# The path to this module
module_base <- file.path(repository_base, "analyses", "cell-type-wilms-tumor-06")In this notebook, we are working with the Wilms tumor sample defined
in SCPCS000194 from the Wilms tumor dataset SCPCP000006. We work with
the pre-processed and labeled Seurat object that is the
output of
02b_label-transfer_fetal_kidney_reference_Stewart.Rmd saved
in the results directory.
result_dir <- file.path(module_base, "results", params$sample_id)
# predictions files
predictions_file <- glue::glue("{params$sample_id}_copykat_prediction.txt")
predictions_paths <- c(
no_ref = file.path(result_dir, "copyKAT", "noref", params$distance, predictions_file),
with_ref = file.path(result_dir, "copyKAT", "ref", params$distance, predictions_file)
)
# full results gene by cell
full_ck_result_file <- glue::glue("{params$sample_id}_copykat_CNA_results.txt")
full_ck_result_paths <- c(
no_ref = file.path(result_dir, "copyKAT", "noref", params$distance, full_ck_result_file),
with_ref = file.path(result_dir, "copyKAT", "ref", params$distance, full_ck_result_file)
)Reports will be saved in the notebook directory. The
pre-processed and annotated Seurat object per samples are
saved in the result folder.
Here we defined function that will be used multiple time all along the notebook.
Seurat objectBelow we look at the heatmaps produced by CopyKAT.
Below we prepare and plot a UMAP that shows which cells are
classified as diploid, aneuploid, and not defined by
CopyKAT. We show a side by side UMAP with results from
running CopyKAT both with and without a reference of normal
cells.
# read in ck predictions from both reference types (no_normal and with_normal)
ck_results_df <- predictions_paths |>
purrr::map(readr::read_tsv) |>
dplyr::bind_rows(.id = "reference_used")
# get umap coordinate
umap_df <- srat[["umap"]]@cell.embeddings |>
as.data.frame() |>
rownames_to_column('barcodes')
cnv_df <- umap_df |>
dplyr::left_join(ck_results_df, by = c("barcodes" = "cell.names")) ggplot(cnv_df, aes(x = umap_1, y = umap_2, color = copykat.pred)) +
geom_point(alpha = 0.5, size = 0.5) +
theme_bw() +
facet_wrap(vars(reference_used))To validate some of these annotations, we can also look at some commonly found copy number variations in Wilms tumor patients:
Although these are the most frequent, there are patients who do not have any of these alterations and patients that only have some of these alterations. See Tirode et al., and Crompton et al..
CopyKAT outputs a matrix that contains the estimated
copy numbers for each gene in each cell. We can read that in and look at
the mean estimated copy numbers for each chromosome across each cell. We
might expect that tumor cells would show an increased estimated copy
number in Chr1q, and/or a loss of Chr1p, Chr11p and Chr16q.
# read in full gene by cell copy number detection results
full_ck_results_df <- full_ck_result_paths |>
purrr::map(readr::read_tsv) |>
dplyr::bind_rows(.id = "reference_used")
# for every cell, calculate the mean detection level across all genes in a given chromosome
full_cnv_df <- full_ck_results_df |>
tidyr::pivot_longer(
cols = -c(
reference_used,
chrom
),
names_to = "barcodes",
values_to = "cnv_detection"
) |>
dplyr::group_by(chrom, barcodes, reference_used) |>
dplyr::summarise(mean_cnv_detection = mean(cnv_detection))
# join with cnv info
cnv_df <- cnv_df |>
dplyr::left_join(full_cnv_df, by = c("barcodes", "reference_used")) |>
dplyr::filter(!is.na(chrom))Let’s look at the distribution of CNV estimation in cells that are
called aneuploid and diploid by CopyKAT.
# create faceted density plots showing estimation of CNV detection across each chr of interest
# colored by aneuploid/diploid estimation
ggplot(cnv_df, aes(x = mean_cnv_detection, color = copykat.pred)) +
geom_density() +
theme_bw() +
facet_grid(
rows = vars(chrom),
cols = vars(reference_used)
)# record the versions of the packages used in this analysis and other environment information
sessionInfo()## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Vienna
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] copykat_1.1.0 ggalluvial_0.12.5 org.Hs.eg.db_3.19.1 AnnotationDbi_1.66.0 IRanges_2.38.1 S4Vectors_0.42.1 Biobase_2.64.0
## [8] BiocGenerics_0.50.0 clusterProfiler_4.12.6 enrichplot_1.24.4 msigdbr_7.5.1 patchwork_1.2.0 lubridate_1.9.3 forcats_1.0.0
## [15] stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1
## [22] tidyverse_2.0.0 SCpubr_2.0.2 sctransform_0.4.1 Seurat_5.1.0 SeuratObject_5.0.2 sp_2.1-4 optparse_1.7.5
## [29] Matrix_1.7-0
##
## loaded via a namespace (and not attached):
## [1] fs_1.6.4 matrixStats_1.3.0 spatstat.sparse_3.1-0 httr_1.4.7 RColorBrewer_1.1-3 tools_4.4.1 utf8_1.2.4
## [8] R6_2.5.1 DT_0.33 lazyeval_0.2.2 uwot_0.2.2 withr_3.0.1 gridExtra_2.3 progressr_0.14.0
## [15] cli_3.6.3 spatstat.explore_3.3-2 fastDummies_1.7.4 scatterpie_0.2.4 labeling_0.4.3 sass_0.4.9 spatstat.data_3.1-2
## [22] ggridges_0.5.6 pbapply_1.7-2 yulab.utils_0.1.7 gson_0.1.0 DOSE_3.30.5 R.utils_2.12.3 parallelly_1.38.0
## [29] rstudioapi_0.16.0 RSQLite_2.3.7 generics_0.1.3 gridGraphics_0.5-1 vroom_1.6.5 ica_1.0-3 spatstat.random_3.3-1
## [36] GO.db_3.19.1 fansi_1.0.6 abind_1.4-5 R.methodsS3_1.8.2 lifecycle_1.0.4 yaml_2.3.10 qvalue_2.36.0
## [43] Rtsne_0.17 grid_4.4.1 blob_1.2.4 promises_1.3.0 crayon_1.5.3 miniUI_0.1.1.1 lattice_0.22-6
## [50] cowplot_1.1.3 KEGGREST_1.44.1 pillar_1.9.0 knitr_1.48 fgsea_1.30.0 future.apply_1.11.2 codetools_0.2-20
## [57] fastmatch_1.1-4 leiden_0.4.3.1 glue_1.7.0 ggfun_0.1.6 spatstat.univar_3.0-0 data.table_1.16.0 vctrs_0.6.5
## [64] png_0.1-8 treeio_1.28.0 spam_2.10-0 gtable_0.3.5 cachem_1.1.0 xfun_0.47 mime_0.12
## [71] tidygraph_1.3.1 survival_3.7-0 fitdistrplus_1.2-1 ROCR_1.0-11 nlme_3.1-166 ggtree_3.12.0 bit64_4.0.5
## [78] RcppAnnoy_0.0.22 GenomeInfoDb_1.40.1 rprojroot_2.0.4 bslib_0.8.0 irlba_2.3.5.1 KernSmooth_2.23-24 colorspace_2.1-1
## [85] DBI_1.2.3 tidyselect_1.2.1 bit_4.0.5 compiler_4.4.1 httr2_1.0.3 plotly_4.10.4 shadowtext_0.1.4
## [92] scales_1.3.0 lmtest_0.9-40 rappdirs_0.3.3 digest_0.6.37 goftest_1.2-3 spatstat.utils_3.1-0 rmarkdown_2.28
## [99] XVector_0.44.0 htmltools_0.5.8.1 pkgconfig_2.0.3 highr_0.11 fastmap_1.2.0 rlang_1.1.4 htmlwidgets_1.6.4
## [106] UCSC.utils_1.0.0 shiny_1.9.1 farver_2.1.2 jquerylib_0.1.4 zoo_1.8-12 jsonlite_1.8.8 BiocParallel_1.38.0
## [113] GOSemSim_2.30.2 R.oo_1.26.0 magrittr_2.0.3 GenomeInfoDbData_1.2.12 ggplotify_0.1.2 dotCall64_1.1-1 munsell_0.5.1
## [120] Rcpp_1.0.13 ape_5.8 babelgene_22.9 viridis_0.6.5 reticulate_1.38.0 stringi_1.8.4 ggraph_2.2.1
## [127] zlibbioc_1.50.0 MASS_7.3-61 plyr_1.8.9 parallel_4.4.1 listenv_0.9.1 ggrepel_0.9.5 deldir_2.0-4
## [134] Biostrings_2.72.1 graphlayouts_1.1.1 splines_4.4.1 tensor_1.5 hms_1.1.3 igraph_2.0.3 spatstat.geom_3.3-2
## [141] RcppHNSW_0.6.0 reshape2_1.4.4 evaluate_0.24.0 renv_1.0.7 BiocManager_1.30.25 tzdb_0.4.0 tweenr_2.0.3
## [148] httpuv_1.6.15 RANN_2.6.2 getopt_1.20.4 polyclip_1.10-7 future_1.34.0 scattermore_1.2 ggforce_0.4.2
## [155] xtable_1.8-4 RSpectra_0.16-2 tidytree_0.4.6 gitcreds_0.1.2 later_1.3.2 viridisLite_0.4.2 aplot_0.2.3
## [162] memoise_2.0.1 cluster_2.1.6 timechange_0.3.0 globals_0.16.3